{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Yijun Wang\n",
    "### 2020/02/09\n",
    "\n",
    "\n",
    "## Note for parameter\n",
    "- Generation period value is 7.5, taken from reference 1 and 2\n",
    "- P value is 0.695, taken from People's Daily Weibo. https://m.weibo.cn/u/2803301701\n",
    "- Formula to estimate R0 refers to reference 3\n",
    "\n",
    "## Reference\n",
    "- 1. Li, Q., Guan, X., et al. (2020, January 29). Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus–Infected Pneumonia. The New England Journal of Medicine. https://www.nejm.org/doi/full/10.1056/NEJMoa2001316#article_references\n",
    "\n",
    "- 2. Wu, Joseph T., et al. (2020, January 28). Nowcasting and Forecasting the potential domestic and International Spread of the 2019-nCoV Outbreak Originating in Wuhan, China: a Modeling Study. Lancet.https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30260-9/fulltext\n",
    "\n",
    "- 3. Zhou Tao, Liu, Y., et al. (2020, January 29). Preliminary Prediction of the Basic Repreduction Number of the Novel Coronavirus 2019-nCoV http://kns.cnki.net/kcms/detail/51.1656.r.20200204.1640.002.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load data\n",
    "import requests\n",
    "import pandas as pd\n",
    "\n",
    "url = 'https://lab.isaaclin.cn/nCoV/api/overall?latest=0'\n",
    "r = requests.request('GET', url)\n",
    "data = r.json()\n",
    "df = pd.DataFrame.from_records(data['results'])\n",
    "from datetime import datetime\n",
    "import pandas\n",
    "\n",
    "# date\n",
    "df['t'] = pandas.to_datetime(df['updateTime']/1000,unit='s')\n",
    "import datetime\n",
    "df = df.resample('D', on = 't').max()\n",
    "df['date'] = df['t'].dt.date"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "# R0\n",
    "import math\n",
    "import numpy\n",
    "from array import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# generation period or serial interval (Wu, 2020; Li, 2020)\n",
    "Tg = 7.5\n",
    "\n",
    "# days: 疾病已爆发时间\n",
    "import datetime\n",
    "df['days'] = (df['updateTime']/1000 - datetime.datetime(2019,12,1,0,0,0).timestamp())/60/60/24\n",
    "\n",
    "def R0Func(confirm, suspect,t):\n",
    "    # confirm是确诊人数；susp是疑似人数；t是疾病已爆发时间\n",
    "    # p为疑似病例转化为确诊病例的概率\n",
    "    p = 41/59 # 人民日报微博\n",
    "    # yt为实际预估感染人数\n",
    "    Yt = suspect * p + confirm\n",
    "    # lamda预估感染人数的增长率\n",
    "    lamda = math.log(Yt)/t\n",
    "    R0 = 1 + lamda * Tg + p * (1 - p) * pow(lamda * Tg, 2)\n",
    "    return R0\n",
    "\n",
    "df['R0'] = df.apply(lambda x: R0Func(x['confirmedCount'], x['suspectedCount'], x['days']), axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAIBCAYAAABusDQ8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUZdrH8e8QEkLoJUIECTUgRUIvui8gKEgTBBVsSAuCqCgCrqAgChbcLKuuioJGdAUUBAUEFlAUVxQsEaUkoUSKQek9oeR5/3iSgYF0MpnMzO9zXedK5pwz59xnMplzz1MdxhiDiIiIiBco4ukARERERHJKiYuIiIh4DSUuIiIi4jWUuIiIiIjXUOIiIiIiXkOJi4iIiHgNJS4iIiLiNZS4iIiLNWvW4HA4OHfuXIGf2+FwsGrVqgI/78XeffddrrnmGooUKUJMTIxHYxGRyylxEfFy7du3x+FwXLbMnTs32+dOmDCB9u3bu6xr27YtSUlJFC1a1E0Rw6pVq3A4HJetT0pK4v/+7//cdt7snDlzhhEjRjBu3Dj27t3LnXfeedk+MTExztc4ICCAatWq8cgjj3Dq1CmX/eLj4+nQoQPFixenevXqvPPOOwV1GSI+zX2fTCJSYEaNGsW4ceNc1pUtWzZPxwoKCqJy5cr5EVaueeq86f744w+Sk5Pp1q0bYWFhme4XFhbGTz/9RGpqKlu2bGHgwIEEBgby8ssvA3D27Fm6detGZGQkGzZs4Pvvv2fYsGGEh4fTsWPHgrocEZ+kEhcRH1CiRAkqV67ssgQHBwOwY8cOunTpQunSpSldujStWrVi27ZtxMTEMGXKFL766itnCUJiYuJlVUWTJk3ihhtu4F//+hdhYWGULVuWqVOnkpKSwgMPPEDp0qWpXbs2K1eudMazdetWunbtSsWKFSlbtixdu3Zl586dACQmJnLTTTcBOM+bXiVzaVXRsmXLaNSoEcWKFaN27drMnj3buS0xMRGHw8GiRYto2bIlJUqUoH379uzatSvL12r27NnUrl2bYsWK0ahRI5YtWwbYKrIaNWoAULNmTefrkZEiRYpQuXJlrr76ajp27Mjtt9/O6tWrXeLevXs377zzDg0bNmTw4MH079+fV199Ndu/pYhkTYmLiI8bOXIklSpVYsOGDfzwww88/PDDFClShDvvvJNRo0bRpk0bkpKSSEpK4pprrsnwGBs3biQ2NpYvv/yS6dOnM378eHr27EmDBg348ccf6dy5M/fddx9nzpwB4MSJE/Tt25dvvvmGb775hqCgIPr16wfANddcw0cffQTgPG9GVTKJiYn06tWLXr16sXHjRkaNGsWgQYP43//+57LfpEmTePHFF1m/fj2nTp3i0UcfzfS1+Pbbbxk0aBAPP/wwGzdupHfv3vTq1YvExETatm3LunXrAFi/fn2Wr8fFdu3axYoVKwgKCnKuW79+PS1atKBUqVLOdR07duT777/P9ngikg0jIl6tXbt2JjAw0JQoUcJl2b59uzHGmIYNG5rZs2dn+Nzx48ebdu3auaz78ssvDWDOnj1rjDFm4sSJply5ciY5Odm5T926dU23bt2cj5OSkgxgNm7cmOF50rf//vvvxhhjVq5caTL6+AHMypUrjTHGjBs3zrRo0cJl+5133mn69u1rjDFm586dBjDz5s1zbv/www9NhQoVMowh/fm33367y7pWrVqZxx9/3BhjTEJCggHMzp07Mz3Gu+++axwOhylRooQpXry4AYzD4TBz58517jN06FBz2223uTxv6dKlJiAgINPjikjOqMRFxAcMHTqU2NhYlyW9tGDEiBEMGTKEzp078/LLL7N79+5cH79OnToUK1bM+bhSpUo0aNDA5THA/v37ATh69CgPPvggderUcVYlAbk6d1xcHK1bt3ZZ16ZNG+Li4lzWNWrUyPl75cqVOXjwIOfPn7+iY2anUqVKxMbGsn79eh599FH69+/vUmpkjMnV8UQk55S4iPiAcuXKUbt2bZclMDAQgOHDh7Nlyxa6du3K8uXLqVevHmvXrs3V8dOPlc7hcLisS+8hlJqaCsDo0aP56quvmD59Ot999x3ffvstYBut5lROb/4ZxZHZc/MroQgICKB27do0bNiQ6Ohodu3axYwZM5zbK1WqxF9//eXynP379xMaGpov5xfxZ0pcRPxAzZo1eeSRR1i1ahXt2rVjzpw5gL3pZ1Y6cSW+++47hgwZQrdu3ahfvz5Hjx512Z6ebGR17nr16vHdd9+5rFu3bh316tXLc1zuOCbAE088wdNPP+3sEt2yZUt++OEHTpw44dzniy++oFWrVld0HhFR4iLiE06ePMm+fftclpMnTwLw6KOPsmrVKhITE1m7di0bN26kbt26AISHhxMXF8fWrVs5cOCAs8TkStWqVYv58+ezefNmvvnmG8aMGeOyPTw8HIDPP/+cAwcOkJKSctkxhg8fzi+//MLTTz9NfHw8r732GvPnz2fUqFF5juvhhx/mk08+4bXXXiM+Pp6nn36an3/+mREjRuT5mAC33HILpUqVcpa6dOnShSpVqjBo0CA2bdrEO++8w5w5c3jooYeu6DwiosRFxCdMnz6dsLAwlyW96+3Zs2eJioqiXr169O/fn7vuuouRI0cC0LdvX1q2bEmLFi0IDQ3NtitxTv3jH//AGEOzZs2Iiopi8uTJLturV6/OuHHjGDhwIKGhoc4SoIuFh4ezaNEiFi5cSMOGDZk+fTqzZs2ibdu2eY6rbdu2vPPOO0yfPp2GDRuycOFCFi1aRPXq1fN8TLDdo0eMGMHLL79MSkoKQUFBLF26lL/++otmzZrxzDPP8MYbb2gMF5F84DBqRSYiIiJeQiUuIiIi4jWUuIiIiIjXUOIiIiIiXkOJi4iIiHgNJS4iIiLiNZS4iIiIiNdQ4iIiIiJeQ4mLiIiIeA0lLiIiIuI1lLiIiIiI11DiIiIiIl5DiYuIiIh4DSUuIiIi4jWUuIiIiIjXUOIiIiIiXkOJi4iIiHgNJS4iIiLiNZS4iIiIiNdQ4iKFxoMPPkj58uVxOBwkJiZ6Opxs3X///dxzzz35drw1a9bgcDg4d+5cvh2zsJgwYQLt27fP1XNuuOEGJk2a5JZ4svLuu+9yzTXXUKRIEWJiYgr8/HkRExND1apVPR2GT7+HpfBQ4iI5FhMTQ0BAAJMnT75s2/3334/D4cDhcFC8eHEaNWrEBx984LI9q5v8119/zdtvv83SpUtJSkrimmuuccs15FT79u2d1xMQEEDVqlV55JFHSElJ8WhceTFp0iRuuOEGT4fhdlWrVr3iROPMmTOMGDGCcePGsXfvXu688878CS6PJk2a5HwfBgUFERYWRvfu3fnss89c9rvzzjv5+eef3RrLtGnTuOaaazh27Jhz3ZkzZ2jQoAGjRo0CoG3btiQlJVG0aFG3xpKZhIQEQkJC+Pzzz13W79u3j/Lly/Pee+95JC7JX0pcJMdmz57NqFGjmD17dobb+/TpQ1JSEps3b6Zv377cd999rF27NkfH3rFjB2FhYbRp04bKlSsTEBBw2T5nzpy5ovhza9SoUSQlJbFr1y5iYmL45JNPePbZZws0hiuRmpqqb7659Mcff5CcnEy3bt0ICwujePHil+1T0Mlry5YtSUpKYufOnSxevJimTZtyxx138MQTTzj3KV68OKGhoW6N47HHHqNy5cqMGTPGue65557j9OnTTJkyBYCgoCAqV67s1jiyUqdOHSZPnsywYcM4fvy4c/3IkSNp1aoVAwYM8Fhskn+UuEiO7Nq1ix9++IFnn30WYwzffPPNZfsEBwdTuXJlatSowcSJE6lduzZLly7N9tiTJk1i4MCB7Nq1C4fDQfXq1QGoXr06L774In369CEkJIRXXnkFgGXLltGoUSOKFStG7dq1XRKpxMREHA4Hn3zyCc2bN6d48eJ06tSJgwcP8vHHH1OrVi3KlSvHo48+ijEmy7hKlChB5cqVqVKlCp06daJPnz5Zfqt9/vnnufbaawkJCaFOnTrOeNO1b9+esWPHMmzYMEqVKkX16tWZO3dupsfbt28fDRs2ZNiwYc5YV61a5byuiIgI/v3vf1927fPnz6dly5YEBwcTGxub5TW66zo+/vhjwsPDKVGiBAMGDCA5OTnLcxpjGD9+POXKlSM0NJRp06Zdts+oUaOoWbMmISEhNGjQgHnz5rnEtHfvXgYOHIjD4XBWS3322We0bt2aUqVKcfXVVzNixAhOnjyZYQxr1qyhRo0aANSsWdNZZXn//fdz991388QTT1CxYkX69u0L2G/3N998M8WLF+eqq65izJgxLoli9erVmTZtGrfddhshISHUr1+fDRs28Ouvv9KqVStKlixJt27dOHToUJavTWBgoPN92Lx5cyZPnszMmTN56aWX2LRpE3B5VVF6CeeECRMoX748V199NdHR0S7HXbVqFXXr1qV48eJ0796dF1980fm/l5GAgADeeecdYmJi+PLLL9m4cSMvvvgib7/9NiVKlHC+hhdXFX377bd06NCBsmXLEhoaSv/+/Tlw4IDzmOlxz58/nxo1alC2bFkGDRrkkhwmJibSvn17goODiYyM5KOPPsqyOvnRRx/l6quvZuzYsQAsXLiQlStX8tZbbzmv2+FwuDxn5syZLte+aNEiWrVqRalSpahSpQojR47k1KlTzu3pVZ//+te/qFy5MqGhoS6JpLiZEcmBZ5991tx1113GGGP+/ve/myFDhrhsHzBggLn77rtd1l133XXmsccey3R7uuPHj5t//OMfpmrVqiYpKcn89ddfxhhjwsPDTfny5c1bb71ltm/fbnbv3m127txpgoKCzIQJE8zWrVvNq6++agICAsw333xjjDFm586dBjCNGzc2a9asMT///LOpXbu26dChg+nZs6f59ddfzZIlS0xQUJD57LPPMr3edu3amfHjxzsf79q1yzRo0MBMnDgx02v+xz/+YdauXWt27Nhh5s2bZ0qUKGGWLl3qcszSpUub6Ohok5CQYCZOnGiCg4PNn3/+aYwx5ssvvzSAOXv2rNm9e7eJiIgwDz/8sElNTTXGGLN161ZTqlQpM3PmTLN9+3azePFiExoaaubOnety7fXq1TMrVqwwCQkJ5siRI2bixInm+uuvz/RaL3Wl17Ft2zZTtGhRM3nyZLN161YzefJkU7JkSdOuXbtMzxkTE2NKlixpPvroI/Pbb7+Z3r17m5IlS7q83pMnTzbff/+92b59u3njjTdMYGCg2bhxozHGmIMHD5qwsDAzffp0k5SUZA4ePGiMMWbevHnms88+M9u3bzdfffWVqVevnhkzZkyGMaSkpJh169YZwKxfv94kJSWZc+fOmQEDBpgSJUqYhx9+2GzdutXEx8ebc+fOmWuvvdZ0797dbNy40Xz++efmqquuMlOmTHEeLzw83Fx11VXm/fffN3FxcaZXr16mbt26pkOHDi7vzfT/kYxk9rc7f/68qVChgnnhhReMMca8++67pkqVKs7tAwYMMKVKlTJjx441cXFxZsaMGQYwv/zyizHGmEOHDplSpUo5r2nGjBmmXLlyJjw8PNNY0j399NOmVq1apmnTppd9Dlz8HjbGmBUrVph58+aZhIQEs2HDBnP99deb22+/3bn/u+++a4KDg52v4xdffGHKly9vXnnlFec+119/vWnXrp2JjY01K1euNBEREQYwO3fuzDTGTZs2meDgYLNo0SITFhZm3nzzTee2lStXmktvfW+//bbLtc+ZM8csWbLEbN++3axZs8ZERESYv//9787t48ePN6VKlTIDBw40W7ZsMQsWLDBFixY1n3/+ebavn1w5JS6SI3Xq1HHe6H/55RdTunRpc+rUKef2i2/i58+fNx988IEBzKeffnrZ9oxc+sFhjP3gv//++13WjRs3zrRo0cJl3Z133mn69u1rjLlw8543b55z+/PPP28cDofzxmqMMZ07d87yhtGuXTsTGBhoSpQoYYKDgw1gbrrpJnPmzJkMrzkjw4YNMwMHDnQ55i233OJ8fPbsWRMSEmIWL15sjLnwoZ+QkGBq1Khx2Q124MCBZvTo0S7rpkyZYjp27Ohy7TExMS775DZxudLrGDt2rGnVqpXLMVq1apVl4tKyZUszbtw45+NDhw6Z4sWLuyQul+rcubN55plnnI+rVKli3n333SyvZc6cOaZGjRqZbk9ISLjspjhgwABTs2ZNc/78eee6ZcuWmeDgYGeCZIwxb7zxhqlYsaLzcXh4uBk+fLjzcXpS9PHHHzvXPf/886Zp06aZxpPV365169bmgQceMMZknLjUr1/fZf+IiAjz6quvGmOMee2110y1atVcrql///45SlxSUlJM5cqVTYUKFcyRI0dctl2auFxq3bp1pmjRoubcuXPOuB0Oh9m3b59zn6ioKNOnTx9jjDG//fab838iXXoSllXiYoz9slWkSBHToUMHZ/JvTM4Sl0u9//77pk6dOs7H48ePNxUrVjQpKSnOdTfeeKPLe1jcR1VFkq1vv/2W/fv307lzZwCuu+46qlatyqeffuqy37x58yhZsiTBwcEMGzaMZ555hp49e17RuZs0aeLyOC4ujtatW7usa9OmDXFxcS7rGjVq5Py9UqVKhIaGctVVV7ms279/f5bnHjp0KLGxsfzyyy8sX76cPXv2OBshZmTp0qXccMMNVKpUiZIlS/LOO++we/fuTOMqWrQoFStW5K+//nLZp3379nTv3p2XXnrJZf2vv/7Ka6+9RsmSJZ3L5MmT2bFjh8t+l75muXWl1xEXF0fLli1d9r/08aUufU65cuWoXbu2yz7vvfcezZs3p2LFipQsWZLVq1dfFtelNm/eTO/evalWrRqlSpVi4MCB2T4nI40bN6ZIkQsfl3FxcdSpU4fy5cs717Vp04YDBw64VP1c+j4EaNCggcu67N6HmTHGXFblcbGGDRu6PK5cubLzb5SQkEBkZKTLNTVv3jxH512xYgWHDx/m6NGjbN68Oct99+zZw7333kvNmjUpVaoUHTt25Ny5c+zbt8+5T2hoqPO1ySjO0qVLu7wXchrnE088QWpqKmPHjs3ydcrIpk2buPXWW53vm6FDh172vomIiCAoKCjDuMW9lLhItmbPns2RI0cICQmhaNGiFC1alC1btlzWQr9bt27Exsayc+dOjh8/ztNPP33F5w4JCXF5bLJpl5IuMDDQ+bvD4XB5nL4uNTU1y2Ok3zwjIiLo3LkzEydOZMaMGZw+ffqyfXfs2MFtt93GjTfeyNKlS/n555+57777OHv2bKZxZRZHly5d+Pzzz0lKSnJZf+LECR577DFiY2Ody2+//cYXX3zhst+lr1lu5Md1ZHdDzUxWz1m7di1Dhw7l3nvvZeXKlcTGxtKpU6fL4rpUz549cTgc/Oc//+GHH37glVdeyVOD5fx6H2a0Lrv3YUZSU1NJSEjIsk2KO/5GR48eZfjw4UyePJkRI0YwdOjQLP8G999/P7///jtvvfUWGzZsYP78+QAuz3FHnICzZ9OlPZzSk7WL/4aXXkP37t0JCgpyvm/++c9/5ul/WdzDM33WxGukpKQwb948YmJiaNasmXP9X3/9xc0330xSUhJhYWEAlCxZ8rJvyfmtXr16l92o161bR7169dx6XrAfgOfPn+fMmTOX9Tb56aefKF68uEtX8Z07d+bpPK+//jqDBg3ipptu4quvvqJChQqA/dYfFxfn1tc4P66jbt26l/Um27BhA8WKFcv0OREREaxfv57evXsDcOTIEbZt2+bc/v3331O/fn0eeeQRwN64t2/f7vJNPTAwkPPnzzsfHzhwgO3btzN//nwiIyMB+Oijj3J1LZmpV68eCQkJHDp0yFnqsm7dOkJDQ11KYdzlP//5D4cPH85ziWZERASLFi0iNTXVeSP/8ccfs33e448/TuXKlRk9ejTJyck0aNCA559/PtMvKd999x0ffPABnTp1Auz7ILdxHj16lO3bt1OrVq0cx5mV9N5X+/btc352/frrr87t+/btIzExkcWLFztLrT788MMrOqfkL5W4SJbSq4PuuusuGjZs6FxuvPFGrr32WpexWgrC8OHD+eWXX3j66aeJj4/ntddeY/78+VlW4eTVyZMn2bdvH0lJSfzvf//jueee44YbbqBMmTKX7VurVi2OHTtGTEwM27Zt47nnnsv1h3S69IHP6tSpQ+fOnZ3jZowZM4YlS5YwYcIENm/ezKZNm4iJieH111/P0bVcXFITGxvrUlyfn9cRFRXFhg0bmDJlCvHx8UyZMoXffvsty+cMHz6cf//738yfP5/NmzczZMgQly7xtWrVIi4ujiVLlhAXF8dDDz10Wfzh4eF8/fXX7Nu3j6NHj1KuXDnKlSvH22+/zY4dO5g3bx4zZszI1bVk5uabb6ZGjRrcf//9/PbbbyxbtoyJEye65X149uxZ9u3bx969e/nxxx95+umnGTp0KE8++WSeE/a77rqLw4cPM3r0aOLj45k1axbLly/PsnTjiy++YPbs2cyaNYuAgABKlCjBjBkzmDp1Klu2bMnwObVq1eL9998nISGB5cuXM3Xq1FzF2aBBA66//nqGDh3Kxo0bWb16tbN3VF5LYiIiIqhUqRITJ05k27ZtzJ49mwULFji3V6xYkTJlyjjfN3PmzGHmzJl5Ope4hxIXydJ7771Ht27dLisWBbj11lsLfECn8PBwFi1axMKFC2nYsCHTp09n1qxZtG3bNt/PNX36dMLCwqhSpQp9+/a9rAvuxZo0acKUKVMYO3YsTZs2JTExkWHDhuX53EWLFmXu3LmUL1+ebt26cerUKZo1a8bKlSv56quvaNasGTfccAPvvvtultUF6WJjY2nSpInL8uabb7rlOmrXrs0HH3zAjBkzaNKkCZs3byYqKirL59x///08+OCDDBkyhP/7v/+jefPmNG7c2Lm9V69ezqqitm3bUqpUKXr06OFyjEmTJvH9999zzTXXcOuttxIQEMB//vMf/vvf/9KgQQNmzJiR4eCJeVGkSBE+/fRTTp8+TYsWLRgwYAD33Xefswtuflq/fj1hYWFUr16dbt268dNPP/Hxxx/z3HPP5fmY5cqVY8GCBSxdupTGjRuzcOFCHn744UxLxU6dOsWQIUMYO3asy9+lc+fO3HHHHQwdOjTD6rOZM2eybds2GjVqxFNPPZWnmN9//33Onz9Py5Yteeyxx5zdjrMqwctKsWLF+OCDD/j666+57rrr+PTTT13+bkWLFuWDDz5g6dKlNGjQgFmzZvHMM8/k6VziHg6T08paERHxWUOGDCEpKSlHYy950gcffMDw4cM5evSoS+Ni8R9q4yIi4odiYmKoV68eoaGhrFy5kvfff79Qzs20cuVKzp07R7169di8eTNPPvkkd911l5IWP6bERUTED+3atYsJEyZw4MABatSowb/+9S/69+/v6bAuk5yczJgxY/j9998JDQ2ld+/evPDCC54OSzxIVUUiIiLiNVTWJiIiIl5DiYuIiIh4DSUuIiIi4jV8unFusWLFnKMkepOUlJQ8j1FQGHhz/N4cO7g3fr02nuPNsYN3x+/NsYP3xr9//35SUlIy3ObTiUtoaCh79uzxdBi5tmLFCueEht7Im+P35tjBvfHrtfEcb44dvDt+b44dvDf+qlWrZrpNVUUiIiLiNZS4iIiIiNfw6aoiERERX5KamprhvFBZuXjW9MLC4XDkefRjJS4iIiKFXGpqKr///jvJycm5el5oaCjx8fFuiurKBAcHEx4enusERomLiIhIIffXX39RpEgR6tSpg8PhyPHzjh07RunSpd0YWd4YY9i7dy9//fUXlStXztVzlbiIiIgUYsYYjhw5QvXq1SlaNHe37SJFihAQEOCmyK5MpUqVSExMpFKlSrlKxtQ4V0REpBAzxmCMITAw0NOh5KvAwEDnteWGEhcREZFCzNfnQlbiIiIi4u+SgWMZLLlr25up6tWrU69ePSIjI6lbty4vvPCCc9uSJUuoV68etWvXpk+fPpw4cSJ/TppGiYuIiIgvSQaqAGWgTLUyUIYLSxXyLXmZP38+sbGxfPnll7zwwgusX7+eEydOMHjwYBYtWsS2bdsICwtjypQp+XPCNEpcRMR3XfStM+BkQL5/6xQplM4AhzLZdihtez66+uqrqVu3Lr///jvLli2jefPm1KtXD4ARI0YwZ86cfD2fehWJiG9K/9aZ9gHeiU4XtpUH9gLBBR+WyBXrCWzPYnt24821BDLraFQL+Cx34WzdupUDBw7Qvn17Zs+eTXh4uHNb9erV2bt3L6mpqXkecO5SSlxEJHPJOL+dOUssAIIo/Df9nHzrLOzXIFKI9e3bF4fDQVxcHP/85z8JDQ0FyFXX5rxQ4iLizTdnd1KJhUjhlF2JyDFse5bMrAfyYUy6+fPn07BhQ1atWkWPHj248cYbqVatGl988YVzn8TERKpUqZJvpS2gNi7i7y5qxEYZ6NSnk1sasXmlAq4nz3eZxS4i+apTp04MHz6cCRMm0KVLFzZs2MDWrVsBeP311+nXr1++nk8lLuLfVJ2QMQPEeTqIPEoE/gm87eE4RDwlCFsqmtFnW/m07fnsqaeeonbt2sTHxzNz5kx69erFuXPnaNSoEe+9916+nkuJi4hYfwCrgVVpyx/Z7D8eeAiIcHNcORULTAPmYRsnNgJ+zWL/R7EJTuGbxkXkygRjq3LPwNGjRylT5qJ6o3yqAk9MTHR5XK5cOQ4ePAhAs2bN6Nmz55WfJBOqKhLxV8eAxcAjQANs1dh9wGygLDAsm+e/BtQFOgEf45mqIwN8AXQGmgAfAh2A/wLfY79dZqQo8A7QEFju/jBFClwwNim/dPGBEmQlLiJZ2ebpAHIgp2OVnAHWAhOBG7A39Z7AK8ARLiQte4FNwEvZnHcJcCfwNXAHUA1bCpN4xVeUvfPAR0ALoCO2hOhO4AdgJXATUBx7LUftsmrBKufvHAXeAA4DtwAD034XkUJPVUUiWfk/4C3gLk8Hkonsev6swCYrq4CvgJNp20oB3bClJZ2AesClPRizqyfvmHaMP4F3gRnAVOB5bDLwANCVzMeLyIvTQAzwMrADm5w8CDwG1LtHM7EAACAASURBVMxg/2Cc3zDPlzjvWi2UHl9U2jGXA28Ct+ZjvCKS71TiIv4tCFstkpFSQEngbmAIcKqggsqF7BoXt8De1FcCzYDJwLdp2z7FtlG5lsuTFrhQT55RicXFXaErAU9gB8Rahi3FWZ72swbwLNm3l8nOwbTjhAMjsCVETwO/Y6usMkpacqIaNuYYbBLYC+gP7L+ycEXyU/q4KL422WL69eR23BeVuIh/CwZGAs8BH8OqIqvo1Cmt1CIIW+VyLzAL22biI+yN3ls8iC0V+Rs2CcutrEosLlUE6JK27MG+Zm9jE4xnsCUZw7AlPOlfmS4aQ8dFegPC37nQQ+gUNnF5BRgElMjD9WTEAQwAbgaGA3OxJVSvYavA3DuWlki2ihQpQmBgIAcPHqRChQq5utGnpqZy/nx2Q+kWPGMMBw8eJDAwMNdjvChxEf9mgDlAZaAXnF99yc05GPuN/AXgKaA58Dr2RudphuwHopqKZ3rNVMW2pRkPLMVWwSwEPsGWjgzDlmxEknGJURlsddPH2PYsjYFxwO2471MrLC3Gj7DJbD9sD6V/p20T8aBq1aqxa9cuDh3K3QBFp0+fpnjx4m6K6soEBgZSrVq1XD9PiYv4t/9hqzgeJ/P/hiLAk9hSi/7A/dieLK+Tf9/6c8Ngq37GYxujFmZFsSUtt2LbpLyNLYkZB0wAzmbyvKPYko+OwFhsY9uCKPlwYBv53oitRpsHrAGmY0veVPoiHhIUFETt2rVJTU3NVZXRqlUXlSIXIg6HI8+j6aqNi/i39HGRclKC8jfsWCFdsb1vmpP1OCHu8C22u29nYCO2gam3qIltuLsHm5S0ymb/Ndgqm5sp+IQhFBvjJ0Ax7PujG7C7gOMQuUSRIkUICAjI8QLkav+CWq5kCgAlLuK/TmG/UTfFjueRExWxY5+8BCRgZ1mdiS0FcadfgO7A9dheQoOAeGz7j8zGKnHTCJlXLAhbqrE0m/2aFEAs2emN7Ro+AFtl2ADby8y32kiKeBUlLuK/FgHHyX17lSLAGGwCcRUwFLgn7Vj5LZ4LbUGWYhuLbsZWt4ST854/knflsb2OPse2vUlvYLwDlzF0XBZ/nuNKxM2UuIj/eg/bBqN/Hp/fBvgZ2+33Q2x349j8CY3d2ISoPrbKoivwE7aEqO4l+140Qqaz54+PjJBZqNyCLX0Zhm3j1BCowIVJOS9e/H2CThE3cmvikpycTK9evYiIiCAyMpIuXbpcNr8BwJo1awgJCSEyMtK5nD592rl91qxZ1KlTh1q1ahEVFcW5c+fcGbb4g73Y9hPdse0Z8qo8tuTmn9gRY1tjR2TNa1XCX9g5dGpjq6DaYkt2llI4qk7yU/oAdxkprNVcpbE9pFZjS9syG9vHG2bPFvFSbi9xiYqKIi4ujtjYWLp3705UVFSG+9WvX5/Y2Fjnkt59a+fOnTz11FN88803bNu2jX379jFr1ix3hy2+7gMglfzp1uwARmF7KIVhB0i7E1tdk1NHsN2ta2J7sKTPofMVdnh+X3RJNZfLUtiruW4E1nk6CBH/5NbEJTg4mK5duzoHy2ndujU7duzI1THmz59P7969qVSpEg6HgwceeIA5c+a4I1zxFwbbZqECtgomv7TAVh3dhh1/pCm2u3JW7SBOAS9iE5bngGuA+WnP64zvd7/15ongsusK/xh2dOJjBRCLiB8p0DYur7zyCj169MhwW1xcHE2bNqVFixa8/vrrzvW7du0iPDzc+bh69ers2rXL7bGKD9sAbMXOP5Tf1RFlsYnHa9huv23IvB1EBWzC8gT2Zh2D7V7dB99PWPzBLOwUAuWxXemfA9ZjB9QTkTxzmAKa/GDq1KksXryY1atXExIS4rLt2LFjGGMoU6YMe/bsoWvXrkyYMIE77riDhx56iGrVqjFmzBgANm3aRI8ePTIsuYmOjiY6Otr5+MiRIyxYsMC9F+YGycnJBAd7w1fOjBX2+K997VqqLanGulfXcayO69fh/Iy9VEIpIp+LJOTPkEz3SSmbwo67drC7y25M0JX/K7rztS/sf9fs5Hf8AScD6NQn84G9vn31W8ptKUeFHytQ/pfyFE22IxyeKXWGg00OcrDZQQ40PUBKaEqGzy9ypgiOszaDTUlJoVixYgCYQENqUGq+XUdB8Ob3jjfHDt4b/+DBg9mzZ0/GG00BmDZtmmnWrJk5fPhwjvafOnWqGTlypDHGmJdeesmMGDHCuW3p0qWmXbt2OTpOlSpVch1rYbB8+XJPh3BFCnX8ycaYcsaYhsaY1Ms353vsu40xZLH8kb+nc+drX6j/rjmQ7/GfNsaUNxn/XcunbU+XYoz50hjzd2NM00v2rW+MGWWMWWaMOZmHY3sBb37veHPsxnhv/Fndv91eVRQdHc2cOXNYuXIlZctmPA1vUlISqan2G8Tx48dZsmQJTZrYLhR9+vRh4cKF/PnnnxhjePPNN+nXr5+7wxZftQQ4jG2UWxDVMdnNE+SJKQMkf+SmcXEQ0B47d9SP2N5jH2Lfh4exDbJvwVYr3QRMI+tZv9VjSfyYW+cq2rNnD6NHj6ZmzZp06NABgGLFivH9998zZMgQevbsSc+ePVmwYAFvvPEGRYsW5dy5c9x+++0MHDgQgJo1a/LMM89w/fXXk5qayo033sjgwYPdGbb4shhsy667PRyH+IaLZs/OlVDs+EH9sY3FfwP+C6wAvsZ21ReRDLk1calatWqmk0HNnDnT+fvIkSMZOXJkpscZOnQoQ4cOzff4xM/8iR22vTOa7VcKDwfQKG0ZDZzGdoW/zZNBiRReGjlX/MeH2B4d+TF2S0554yBr4lnFsbNii0iG3FriIlKovIfthnxrAZ4zvR1ERm0SgvCO8Uqk8PmD7NtPifgolbiIf4jFzrDcj4JPFrx5kDXxjKxK6sBO8vhrAcUiUsgocRH/8F7az4KsJhLJq6xm/Z6F7ZV0A2rEK35JiYv4vrPAf4A62EkQRbxBZrN+D8L2PnJgu1DHeCxCEY9Q4iK+bzmwn4Ibu0XE3TpgJ/W8GhgITCLvM5KLeBklLuL73sMmLPd6OhCRfNQA+A47mecz2ARGA9OJH1DiIr7tIPAZcCNQzcOxiOS3MOAr7Czn76X9POrRiETcTomL+La52DYuapQrvqok8CnwALAa22h3t0cjEnErJS7i297DfrBrFFLxZUWB14EXsdMHtMYOASDig5S4iO/aAmwA+qLJDMX3OYCx2FLGA8DfsA3TRXyMEhfxXeljt9zvySBECtid2PFdAoHuwNueDUckvylxEd90HngfqI795iniT/4GrMM2SI8CxqPu0uIzlLiIb1qFnc/lPvQuF/9UF5u8tACmAvcAKR6NSCRf6CNdfFN6NdF9Ho1CxLMqAV9iJxb9EOgMHPZoRCJXTImL+J6jwEJst9BaHo5FxNNKAAuAh7BjvlwPJHoyIJErU9TTAYjku4+BZNQoVyRdAPAvoAYwGttdegHQKIN9g9DM5VKoKXER3xMDFAdu93AcIoWJA3gU22D3bmyJZEbKY2emVvIihZSqisS3bMNOPtcbO5OuiLjqAyzJYvshNOeRFGpKXMS3zE77qSH+RTLX0tMBiOSdEhfxHanYxKUK0NHDsYiIiFsocRHf8RXwO3AvtjGiiOTNfOwXAZFCSImL+I70sVtUTSRyZQZjex6t8XAcIhlQ4iK+4QT2W2JLoJ6HYxEp7IKwvYcyUg4YiZ1dugPQA9hUQHGJ5IASF/ENnwAnUWmLSE4EY7s8H81g+QN4FdgK9MP2QLoOGJq2TcTDlLiIb3gP+y2yn6cDEfESwdghAy5d0sdvqQnMAdZjJ22cCdQGngKOFXSwIhcocRHv9zvwBdCTzIu/RSRvWmDnO1qMnW39OWwC8zpw1nNhif9S4iLe7/20n6omEnEPB9Ad2Ai8jR1z/UGgAbaa1nguNPE/SlzEuxns2C1XYWe+FRH3KQoMARKAZ4Ek7Ei812NHrBYpAEpcxLutw36I3g0EejgWEX9RApgAbAdGYNvB3ADcBsSl7ZOMbQtzDAJOBjh/J7ngwxXfosRFvJvGbhHxnKuAf2O7S98GLMRWHw0DwoAydunUp5Pzd6qg5EWuiBIX8V6ngblAJNDYw7GI+LO6wAJsdVFL4C3gSCb7ahJHuUJKXMR7fYotelZpi0jh0BabvLyf3Y4ieafERbzXe9jGgnd5OhARcXJghyYQcRMlLuKd/gD+C9yCrWcXERG/oMRFvNN/sLPXqppIxPu8imafljxT4iLexwAx2FFyu3s2FBHJQFaTOBbBdqW+BdhXYBGJD1HiIt7nR2Az0B8o5uFYRORyl0ziuGrBqguTOP6FHcTuv9jJGz/3WJTipZS4iPfR2C0ihd9FkzieL3H+wiSOFbDTBnyMneuoGzAKSPFUoOJtlLiId0kBPgSuBZp7OBYRybu+wC/YEXf/BbQCtng0IvESSlyk8Lto6HA+xg5g1Q99QxPxdtWwM09PAn4FmmFLYzRpo2RBiYsUbsnYIcLThwu/N239RDR0uIgvKIr9f/4KCAWigNuxX1BEMqDERQq3M2T+Aaahw0V8xw1ALDZpWYCdxuNrj0YkhZQSFxERKRzKAfOAWdgvJh2Ap4FzngxKChslLiIiUng4gEHAT9hSl2eBdkCiB2OSQkWJi4iIFD51gXXAY8C32Fng53k0IikklLiIiEjhVAz4B7As7fd+wEDgABd6Gl68qLG+X1DiIoVbELbeOyPl07aLiG/rAmxM+xkDVOJCT8OLF/U09AtKXKRwCwbGpv3+NheGDT+KHVI82ENxiUjBqgQsBaaS+QSN6mnoF5S4SOFmgHexw4Tfy4Vhw0ujpEXE3xQBHvR0EOJpSlykcPsKiAfuRxMqioiIEhcp5Gak/Rzq0ShExFvs9nQA4m5uTVySk5Pp1asXERERREZG0qVLFxITE7Pcv379+jRvfmH2vDVr1hASEkJkZKRzOX36tDvDlsJiP3YEzfbYrpEiItlpjf3Co/mOfJbbS1yioqKIi4sjNjaW7t27ExUVlem+48ePp02bNpetr1+/PrGxsc6lePHi7gxZCov3sNPeD/N0ICJSaARhexRmpBS2d9EDQCdgZ0EFJQXJrYlLcHAwXbt2xeFwANC6dWt27NiR4b5r164lISGBe++9N8Pt4mcM8BZQEejt4VhEpPAIxvYoPJrB8hewGTtR4xdAI+A1Mu+FJF6pQNu4vPLKK/To0eOy9SdPnmTUqFG88cYbGT4vLi6Opk2b0qJFC15//XV3hymFwZdAAmqUKyKXC8a1h+HFPQ1LY6uKVmFnm34IW92c4IlAxR0cxpgCqQmcOnUqixcvZvXq1YSEhLhsGzFiBM2bN2fQoEGsWbOGxx9/nB9++AGAY8eOYYyhTJky7Nmzh65duzJhwgTuuOOOy84RHR1NdHS08/GRI0dYsGCBey/MDZKTkwkO9t6+vvkR/3XPX0fYV2GsnbWWU1VO5VNk2dNr75ljFwRvjt+bYwfPxR9wOoCIdyKotrga54udJ+G+BH7v9TsE5PwYeu09Y/DgwezZsyfjjaYATJs2zTRr1swcPnw4w+2NGjUy4eHhJjw83FSqVMkEBQWZ+vXrZ7jv1KlTzciRI3N03ipVquQ5Zk9avny5p0O4Ilcc/5/GmEBjzI35EU3u+P1r76FjFwRvjt+bYzemEMS/xhhTyxiDMaaNMWZLzp/q8divkLfGn9X92+1VRdHR0cyZM4eVK1dStmzZDPfZuHEjiYmJJCYmMnfuXBo1asSmTZsASEpKIjXVVlAeP36cJUuW0KRJE3eHLZ4Ug22Um3k7bhGRnGuHnTLgUeA77ISNLwLnPBmU5JVbE5c9e/YwevRojhw5QocOHYiMjKRVq1YADBkyhM8++yzbYyxYsIBGjRrRuHFjWrduzU033cTAgQPdGbZ4Uiq2UW4oapQrIvknBIgGvgGqA08AbYHfPBiT5ElRdx68atWqmEya0MycOTPD9e3bt3e2bwEYOXIkI0eOdEt8Ugh9CWzHzk+kCRRFJL+1BX4GJgEvA02Bp4FxQKDnwpKc08i5UrhopFwRcbfi2KqidUAd4CmgJRDryaAkp5S4SOHxJ7AQ6AjU9nAsIuL7WgI/AeOBX4EW2NKXYxeWgJMBFx4neypQuZgSFyk83sU2ltNIuSJSUIoBzwHrgfrAs9iRecvYpVOfTs7fqYKSl0JAiYsUDqnA28BVwK0ejkVE/E9TYAO20e75TPY5BJwpsIgkE0pcpHBYDewABqFGuSLiGUHA3z0dhGRHiYsUDumNcod4NAoRESnklLiI5+0DPgVuAmp5OBYRESnU3DqOixQiyWRcNxuEnZjMk9QoV0REckiJiz9IxraGP5TBtvLYKeI9lbykN8qtBPT0UAwiIumCsJ+LGX1eFkX1FIWA/gT+4AwZ/xOC51vJrwR2YhvlatRKEfG0YOyXuaN2WbVglf19OLZk+J+eDE5AiYt42luAA42UKyKFRzBQ2i7nS5y3v0cD1wHPYEfcFY9R4iKek4RtlHszUMPDsYiIZCUY+BBbMnw3diRd8QglLuI572AHeorydCAiIjnQAJiGrd5+yMOx+DElLuIZ6Y1yKwM9PByLiEhOPQh0BWYDcz0ci59S4uIP0lvJZ6Q4nhmp9r/A78Bg1ChXRLyHA1tafBXwALDLs+H4IyUu/iC9lfxTaY9XYwd9awycBr7wQEwzsB8AGilXRLxNJez4U0eBe8h8biNxCyUu/iKYC98MmmD/8RZhS2Luxs4TVFD+ABYDnYHqBXheEZH80hUYCawFXvRwLH5GiYs/iQdCgXJpj6tjW8kfBW4DThVQHOmNcjVSroh4s5ewDXYnAus9HIsfUeLiT+KBiEvWdQYmA79g62uNm2M4j22UGwZ0d/O5RETcqTj2y18AtuT6hGfD8RdKXPzFwbTl0sQF4ElsEvE+8Iab41iBrbIajCacEBHvdx22qmgb8IiHY/ETSlz8RULaz4wSlyLYpKUWMAr3jgqZPlKuGuWKiK94CFt6/Q4w38Ox+AElLv4iPu1nRokLQFngE2wpSF/gTzfEsBdYAtwChLvh+CIinlAE28uoInZAzd2eDcfXKXHxF9klLmCLPN/G9vrph51QLD/NQo1yRcQ3hWFLXA4D96Eu0m6kxMVfxGOraGpls9/d2GLPNcDf8/H854GZQBVsN0IREV/TAzuL9BrgH54NxZcpcfEX8UA1bCv47LwMtE37+XE+nX85tvhUjXJFxJe9DFwLjAd+9HAsPkqJiz9IxTbOzaqa6GJB2ISlEjAQ2JIPMczAvtvUKFdEfFkItou0A7gLOOnZcHyREhd/8Ad2cLk6uXjO1cBHQDLQmyubwn03sBRbRXTNFRxHRMQbRALPY0u6H/VwLD5IiYs/yEnD3Iz8H3YK9zhgEHkfnO4dbKlPVB6fLyLibR4FOmE7PCz0cCw+RomLP8hqDJfsjALuBBZg625z6xy2UW5VbDdoERF/UAR4D6iArSL/w7Ph+BIlLv4gryUuYOtpZwL1gSfI/UzSy4A92H9cNcoVEX9yNfbz8xC2i3SqZ8PxFUpc/EE8EEjeB30riR2crgR2fJfcDK6U3ih3cB7PLSLizXphq8lXA//0cCw+QomLP4jHjt9yJSUedbHFnvuB24GUHDxnF7bEpRu2qkhExB9FYz9D/w7EejgWH6DExdedBXaQt2qiS/XGVhd9j237kp1Z2KJRjZQrIv6sBLaLNEB/bC9PyTMlLr4uEdtANj8SF4BngY7Am0BMFvulN8q9BuiST+cWEfFWTYHngK3A4x6OxcspcfF1V9IwNyNFgTnYhGQ48HMm+y3FtqIfAgTk07lFRLzZ40AH4A1gHnZ8rIuXZM+F5k2UuPi6/E5cAEKx3aNTgduwLeYv9RY2YVGjXBERqwj2s9GB7ehQ5pKlCkpeckCJi69zR+IC0AJ4DVsVdTcuM6EG/xlsG+V2x/4jioiIdRWZD+Z5CDhTgLF4KSUuvi4e2525shuOPQQ7ou5yYPKF1VWXV7X/mBopV0RE8pkSF18Xjy1tcbjh2A5sqUtTbOLyEXAIqiyvYrs/t0HFniIikq+UuPiyk9hRa/O7muhixbkwE+qdQAUIPhxsz1se1dmKiEi+UuLiy7al/XRn4gIQhupsRUTyw/nsd/F3Slx8WXrD3DoejUJERNIFYUujMzOInI1M7seUuPgyd/UoEhGRvAkG9gJHL1n2A32ARdhpUo57KsDCT4mLL0tI+6kSFxGRwiMYKH3JUhE7KN0D2AkZOwIHPBVg4abExZfFYweLK+fpQEREJFsBwOvAU8AG4G/Abo9GVCgpcfFl6V2h3S2rOtvyadtFRCR7DuzwEtOx8xpdn/ZTnJS4+KqDaUtBJC6X1NmuWrDqQr3t3rTtIiKSc48As7Fzvv0N+MGz4RQmSlx8VXr7loJqmHtRne35Eucv1NsqaRERyZt7sY11T2AnZ/zSs+EUFkpcfJV6FImIeL/uwH+xd+suwCeeDacwUOLiq5S4iIj4hr8BX2E7WtwOzPJsOJ6mxMVXxWMbedXydCAiInLFIoFvgGrYCW5f8mw4nuTWxCU5OZlevXoRERFBZGQkXbp0ITExMcv969evT/PmzV3Wz5o1izp16lCrVi2ioqI4d+6cO8P2DfHYN3hxTwciIiL5ojbwP6AhMA4YS+bTrfgwt5e4REVFERcXR2xsLN27dycqKirTfcePH0+bNm1c1u3cuZOnnnqKb775hm3btrFv3z5mzfLzcrLspGIb56qaSETEt1wNfA20BaZhS1/87Lu8WxOX4OBgunbtisPhAKB169bs2LEjw33Xrl1LQkIC9957r8v6+fPn07t3bypVqoTD4eCBBx5gzpw57gzb+/0BnEKJi4iILyqHbbDbBXgHuANI9mhEBapA27i88sor9OjR47L1J0+eZNSoUbzxxhuXbdu1axfh4eHOx9WrV2fXrl1ujdPrqWGuiIhvKwF8CvQHFgJdgWMejajAFC2oE02dOpWEhATefPPNy7aNGTOGBx98kCpVqpCQkHDZ9vQSGwBjMq/Qi46OJjo62vn4yJEjrFix4gojL3jJyclXFHfVpVVpQAN+PP4jB1YU/GQXVxq/J3lz7ODe+PXaeI43xw7eHX+hj/1euPbEtVRbXI2jzY/y47M/crbsWefmQh9/XpgCMG3aNNOsWTNz+PDhDLc3atTIhIeHm/DwcFOpUiUTFBRk6tevb4wx5qWXXjIjRoxw7rt06VLTrl27HJ23SpUqVxy7JyxfvvzKDvCoMQZjzLb8iCb3rjh+D/Lm2I1xb/x6bTzHm2M3xrvj94rYU40xE4393I8wxsQZY47aZeWClc7fzWmPRZhrWd2/3V5VFB0dzZw5c1i5ciVly5bNcJ+NGzeSmJhIYmIic+fOpVGjRmzatAmAPn36sHDhQv7880+MMbz55pv069fP3WF7twQgEAjPbkcREfF6DmAS8Cq2qUA9oIxdOvXp5PydKvhEWxi3Ji579uxh9OjRHDlyhA4dOhAZGUmrVq0AGDJkCJ999lm2x6hZsybPPPMM119/PbVq1eKqq65i8ODB7gzb+8Vjx28psIpAERHxuJHA22TeRfoQcKbgwnEXt97aqlatmmmblJkzZ2a4vn379vzwg+tsUkOHDmXo0KH5Hp9POgvswDbUEhER/3IH4OO3S42c62sSsX361aNIRER8kBIXX6Ou0CIikplDng7gyilx8TVKXEREJDNNgRnAeU8HkndKXHyNEhcREf8VBJTPZFsJoBjwANAS+LaggspfSlx8TTxQEqjs6UBERKTABQN7gaN2WbVglfN3DmDvEWOBX4HrgXux08R4ESUuviYeW9riyG5HERHxScFAabucL3He+TvBQCngRWzi0gX4AKgLvITXdJVW4uJLTgJ7UDWRiIhkrS7wOfAZcBUwDmgELPNkUDmjxMWXbEv7qcRFRESy4wB6AJuAKdgvvl2BnsB2D8aVDSUuvkQNc0VEJLeCgSeBOKAfsBioD4zHluQXMkpcfIkSFxERyauqwBxgDbYqaSp23qO5ZD6NgAcocfEl6YlLHY9GISIi3qwd8BPwGrbEpT/QHtiInaTxWAZLAU7eqMTFl8QDoUDGk3CLiIjkTFHgQex9ZRiwFojE3l/KZLAU4MzTSlx8SQKqJhIRkfxTEXgT+AFoAaRksl8BzjytxMVXHExblLiIiEh+awr819NBWEpcfEVC2k8lLiIi4g6FZGBTJS6+Qj2KRETEDyhx8RVKXERExA8ocfEV8dhivFqeDkRERHxSVjNPl0/bXgCyTVxOnTrFk08+SY0aNShWrBjFihWjZs2aPPnkk5w4caIgYpSciAeqAcU9HYiIiPikS2aedln2pm0vANkmLgMGDOD48eN8/vnnHDlyhCNHjrBkyRKOHz/OgAEDCiJGyU4q6gotIiLud9HM0y5LASUtYIeYydIvv/zCxx9/7LKufv36vPrqq0RE6E5ZKPwBnEKJi4iI+LxsS1wCAgJISEi4bH18fDwBAQFuCUpySQ1zRUTET2Rb4jJt2jT+9re/0aJFC8LDw3E4HOzcuZMffviBmTNnFkSMkh0lLiIi4ieyTVy6d+/O9u3bWbZsGbt27QKgXbt2zJ07N8OSGPEAJS4iIuInsk1cALZs2YLD4WDAgAFUqFCBTZs2cc899/C///2P/fv3uztGyU48EAiEezoQERER98q2jcuLL77ITTfdxLRp02jdujWvvvoqLVq0oHbt2ipxKSzigdqAmhyJiIiPy7bEJSYmhs2bNxMWFsbWrVtp2LAhK1asoGPHjgURn2TnLLAD6ObpQERERNwv2xKX4OBgwsLCAKhXrx4RERFKWgqTROA8at8iIiJ+IdsSl5SUFLZs2YIxxrnu4sf169d3X3SSPTXMFRERP5Jt4nLq1Cm6du3qsi79scPhYMeOHe6JTHJGiYuIiPiRbBOXpCgoVgAAIABJREFUxMTEAghD8kyJi4iI+BHNDu3t4oGSQGVPByIiIuJ+Sly8XTy2tMXh6UBERETcT4mLNzsJ7EHVRCIi4jeUuHizbWk/lbiIiIifUOLizdQwV0RE/IwSF2+mxEVERPyMEhdvlp641PFoFCIiIgVGiYs3iweuAsp6OhAREZGCocTFm6V3hRYREfETSly81UHgEEpcRETEryhx8VZq3yIiIn5IiYu3Skj7qRIXERHxI0pcvJW6QouIiB9S4uKt4rHzE9XydCAiIiIFR4mLt4oHqgHFPR2IiIhIwVHi4o1SsW1cVE0kIiJ+RomLN/oDOIUSFxER8TtKXLyRGuaKiIifUuLijZS4iIiIn1Li4o2UuIiIiJ9ya+KSnJxMr169iIiIIDIyki5dupCYmHjZfuvWrSMyMpLIyEgaNGjAsGHDSElJAWDNmjWEhIQ4t0dGRnL69Gl3hl34xQOBQLinAxERESlYbi9xiYqKIi4ujtjYWLp3705UVNRl+zRu3JgNGzYQGxvLr7/+yv79+5kxY4Zze/369YmNjXUuxYv7eR/geKA2EODpQERERAqWWxOX4OBgunbtisPhAKB169bs2LHjsv1CQkIIDAwE4MyZM5w+fZoiRVSLlaGzwA5UTSQiIn6pQLODV155hR49emS4LTExkcjISCpWrEjp0qVdSmbi4uJo2rQpLVq04PXXXy+ocAunncB5lLiIiIhfchhjTEGcaOrUqSxevJjVq1cTEhKS6X4nTpzgnnvuoV+/fvTr149jx45hjKFMmTLs2bOHrl27MmHCBO64447LnhsdHU10dLTz8ZEjR1iwYIFbrsedkpOTCQ4OznBb6PehNJ3YlN8e+Y29t+wt4MhyJqv4Cztvjh3cG79eG8/x5tjBu+P35tjBe+MfPHgwe/bsyXijKQDTpk0zzZo1M4cPH87R/nPmzDHdu3fPcNvUqVPNyJEjc3ScKlWq5DjGwmT58uWZb4w2xmCMWVNQ0eRelvEXct4cuzHujV+vjed4c+zGeHf83hy7Md4bf1b3b7dXFUVHRzNnzhxWrlxJ2bJlM9xn+/btnD17FrBtXD755BOuu+46AJKSkkhNTQXg+PHjLFmyhCZNmrg77MJLXaFFRMSPuTVx2bNnD6NHj+bIkSN06NCByMhIWrVqBcCQIUP47LPPANvluUmTJjRu3JgmTZpQqVIlnnrqKQAWLFhAo0aNaNy4Ma1bt+amm25i4MCB7gy7cIsHSgKVPR2IiIhIwSvqzoNXrVoVk0kTmpkzZzp/Hzx4MIMHD85wv5EjRzJy5Ei3xOeV4rGlLQ5PByIiIlLw1OfYm5wE9qBqIhER8VtKXLzJtrSfSlxERMRPKXHxJmqYKyIifk6JizdR4iIiIn5OiYs3SU9c6ng0ChEREY9R4uJN4oGrgIyHwxGR/2/v3oOjqu//jz/DNZKI8RosTAh3E8BkEykiWkvRghAQsf6wVUcNmCqW2irSodQpYGWoKP2hFYuDWIsl6iARwUjH0qJfECtYQUvBoHIxKLcEws0ECJ/vHwn5EhPCQnL2s+/weswwIXuWw5PjZPftnpuINHoaXCw5fiq0iIjIWUqDixVFQDEaXERE5KymwcUKHZgrIiKiwcUMDS4iIiIaXMzQ4CIiIqLBxYyNVNyfqJPvEBEREX80uFhRACQBsb5DRERE/NHgYsExKj5x0W4iERE5y2lwseAr4BAaXERE5KynwcUCHZgrIiICaHCxQYOLiIgIoMHFBg0uIiIigAYXGwqA5kB73yEiIiJ+aXCxoADoDDT1HSIiIuKXBpdodwT4Au0mEhERQYNL9NsElKPBRUREBA0u0U8H5oqIiFTR4BLtNLiIiIhU0eAS7TS4iIiIVNHgEu0KgHOBRN8hIiIi/mlwiXYbgS5AjO8QERER/zS4RLODQCHaTSQiIlJJg0s0+6zyqwYXERERQINLdNOBuSIiItVocIlmGlxERESq0eASzY4PLl28VoiIiEQNDS7RrAC4BEjwHSIiIhIdNLhEswK0m0hEROQEGlyiVRFQjAYXERGRE2hwiVY6MFdERKQGDS7RSoOLiIhIDRpcopUGFxERkRo0uESrAiruT9TJd4iIiEj00OASrQqA9kCs7xAREZHoocElGh2j4q7Q2k0kIiJSjQaXKNSyqCV8g66YKyIi8i0aXKJQ3La4it/oExcREZFqNLhEoVbbWlX8RoOLiIhINRpcolBcoT5xERERqY0GlyjUalsraE7FWUUiIiJSRYNLFIrbFgedgaa+S0RERKKLBpdocwTO+foc7SYSERGpRTPfAVKpFDgMfAZNjjWp2E20D2iBLkInIiJSSZ+4RINSoC1wHpBZ+dhTld+3rVwuIiIiGlyiwmGg+CTLiiuXi4iIiAYXERERsSPQwaW0tJRhw4bRtWtX0tPTGThwIJs3b67xvJUrV5Kenk56ejrdu3fnpz/9KWVlZVXLn3/+ebp06UKnTp3Iycnh6NGjQWaLiIhIlAr8E5ecnBw+/fRT1qxZQ1ZWFjk5OTWek5aWxqpVq1izZg2ffPIJu3btYtasWQBs2rSJRx55hOXLl/PZZ5+xfft2nn/++aCzRUREJAoFOrjExsYyaNAgYmJiALjyyiv54osvajyvVatWNG/eHIDDhw/zzTff0KRJRdr8+fO56aabSExMJCYmhnvvvZfc3Nwgs0VERCRKRfQYl6eeeoohQ4bUumzz5s2kp6dz0UUX0bp166pPZrZu3Ur79v93Cdnk5GS2bt0akd6IaQFccJJlF1QuFxEREWKccy4Sf9GUKVNYtGgRS5cupVWrVid93oEDB7j99tu59dZbufXWWxkzZgxJSUk8/PDDAKxbt44hQ4bU+snN9OnTmT59etX3e/fu5bXXXmv4f0wAEj5OoPe43mwetpl1/28dLVu2BMA1dxxrccxz3ekpLS0lNtbmxWcst0Ow/do2/lhuB9v9ltvBbv/IkSMpLCysfaGLgGnTprnMzEy3Z8+esJ6fm5vrsrKynHPOPf7442706NFVy95880137bXXhrWetm3bnnarNw8753DOve/ckiVLfNfUi+V+y+3OBduvbeOP5XbnbPdbbnfObn9d79+B7yqaPn06ubm5vP322yQkJNT6nM8//5wjR44AFce4LFiwgMsvvxyAm2++mby8PHbs2IFzjj/96U/ceuutQWdHlgMWUHGxuV6eW0RERKJYoINLYWEhDz30EHv37qVfv36kp6fTu3dvAEaNGsUbb7wBwLJlywiFQqSlpREKhUhMTOSRRx4BoGPHjkyaNIm+ffvSqVMnLrnkEkaOHBlkduT9B/gcuAldWUdERKQOgd6rqF27driTHEIze/bsqt+PHDmyzmHknnvu4Z577mnwvqixoPLrcK8VIiIiUU//fx8NFgAXAtf4DhEREYluGlx8+wz4GBiK7tUtIiJyChpcfMur/KrdRCIiIqekwcW3PCAeuM53iIiISPTT4OLTV8BKYBBg7/pAIiIiEafBxafXK79qN5GIiEhYNLj4tICK+xAN8h0iIiJigwYXX4qBZcAPgXP9poiIiFihwcWXRUA52k0kIiJyGjS4+LKAiq0/xHeIiIiIHRpcfDgA/A24FrjIc4uIiIghGlx8eAsoQ7uJRERETpMGFx+OXy13mNcKERERczS4RFoZsBjoDbTz3CIiImKMBpdIWwrsB27yHSIiImKPBpdIW1D5VYOLiIjIadPgEklHgYVAD6Cr5xYRERGDNLhE0gpgNzqbSERE5AxpcImk47uJNLiIiIicEQ0ukeKoGFw6AJd7bhERETFKg0ukrAYKqfi0JcZzi4iIiFEaXCJFu4lERETqTYNLJBzfTdQGuNJzi4iIiGEaXCJhPVBAxSX+tcVFRETOmN5GI0G7iURERBqEBpdIWAAkAN/33CEiImKcBpegbQI+AoYCzT23iIiIGKfBJWivV37VbiIREZF60+AStAVAK+CHvkNERETs0+ASpO1U3J/oBuAczy0iIiKNgAaXIC2k4hou2k0kIiLSIDS4BCmPigNyB/sOERERaRw0uARlL7AUuA44z3OLiIhII6HBJSiLgaPATb5DREREGg8NLkFZQMVdoG/0HSIiItJ4aHAJwiFgCXANcInnFhERkUZEg0sQ/gZ8g84mEhERaWAaXIJw/KaKOr5FRESkQWlwaWiHgUVAJpDkuUVERKSR0eDS0P4JlKDdRCIiIgHQ4NLQju8m0uAiIiLS4DS4NKRyKi7znwJc5rlFRESkEdLg0pBWAjvQQbkiIiIB0eDSkLSbSEREJFAaXBqKo2JwSQIyPLeIiIg0UhpcGspHwBYqPm2J8dwiIiLSSGlwaSh5lV+1m0hERCQwGlwaygLgYuAq3yEiIiKNlwaXhrAB+C8wDGjquUVERKQR0+DSELSbSEREJCI0uDSEBUBr4Ae+Q0RERBq3QAeX0tJShg0bRteuXUlPT2fgwIFs3ry5xvP+8Y9/0Lt3b1JTU+nRowcTJkzAOQfA5s2badasGenp6VW/Pv/88yCzT89WYDUwBGjhuUVERKSRC/wTl5ycHD799FPWrFlDVlYWOTk5NZ5z/vnnk5uby3//+19Wr17NO++8Q25ubtXyhIQE1qxZU/WrU6dOQWeH7/XKr7paroiISOACHVxiY2MZNGgQMTEVFza58sor+eKLL2o8LxQK0bFjx6o/k56eXuvzotICIBYY6DtERESk8YvoMS5PPfUUQ4YMqfM527dvZ/78+QwaNKjqsX379tGrVy8yMjKYPHky5eXlQaeGZxfwP1QMLXGeW0RERM4CMe74wSQBmzJlCosWLWLp0qW0atWq1ufs27eP/v378+Mf/5gHH3wQgLKyMkpKSrjkkksoLi5mxIgRXH/99YwbN67Gn58+fTrTp0+v+n7v3r289tprwfyDgLZL2tLj//fg47Ef8/V1XzfYektLS4mNjW2w9UWa5X7L7RBsv7aNP5bbwXa/5Xaw2z9y5EgKCwtrX+giYNq0aS4zM9Pt2bPnpM/Zt2+f69Onj5s8eXKd65o3b57LysoK6+9t27btaXWetkHOuWbOueKGXe2SJUsadoURZrnfcrtzwfZr2/hjud052/2W252z21/X+3fgu4qmT59Obm4ub7/9NgkJCbU+58CBAwwcOJABAwbwyCOPVFu2c+dOjhw5AlR8+rJgwQJCoVDQ2adWAvwd6Aec77lFRETkLBHo4FJYWMhDDz3E3r176devH+np6fTu3RuAUaNG8cYbbwAwY8YMPvjgA/Ly8qpOeX7ssccAWL58OaFQiLS0NDIyMmjTpg0TJkwIMjs8+cBhdNE5ERGRCGoW5MrbtWtXdT2Wb5s9e3bV7ydMmHDSYWT48OEMHx6F08ECKu4CfaPvEBERkbOHrpx7Jr4B3qLihoqXem4RERE5iwT6iUujU0rF7qF84CBwA7CPiivm2jtoW0RExBx94hKuUqAtcB7w48rHflP5fdvK5SIiIhIoDS7hOgwUn2RZceVyERERCZQGFxERETFDg4uIiIiYocFFREREzNDgIiIiImZocAlXC+CCkyy7oHK5iIiIBErXcQlXLLCN2s8e0nVcREREIkKDy+mIRQOKiIiIR9pVJCIiImZocBEREREzNLiIiIiIGRpcRERExAwNLiIiImKGBhcRERExQ4OLiIiImKHBRURERMzQ4CIiIiJmaHARERERM2Kcc853RFBatmzJxRdf7DvjtB04cID4+HjfGWfMcr/ldgi2X9vGH8vtYLvfcjvY7d+1axdlZWW1LmvUg4tV7dq1o7Cw0HfGGbPcb7kdgu3XtvHHcjvY7rfcDvb7a6NdRSIiImKGBhcRERExo+nEiRMn+o6Qmvr06eM7oV4s91tuh2D7tW38sdwOtvstt4P9/m/TMS4iIiJihnYViYiIiBkaXERERMQMDS4iIiJihgYXERERMUODi0gjtHbtWt8JIiKB0OASZb788kvy8vL47LPPfKeE5Y9//CO7du3ynXHGPvzwQ37xi18wbNgwbrnlFn7729+yY8cO31n1NmTIkMDWvXjx4sDW3ZD+85//sG7dOgA2btzIH/7wB5YuXeq56vTpNSGyGutrgpWf23BocPHsjjvuqPr9smXLuOKKK5g1axbXXHMNeXl5HsvC8/DDD5OcnMzw4cPJz8/H0tn1M2bMIDs7m2PHjrF+/Xouvvhidu7cSSgUYsWKFb7zTmnmzJm1/nrmmWc4ePBgYH/v6NGjA1t3Q3n66acZPHgwAwYMYMaMGYwYMYKCggLuu+8+Zs2a5TuvTnpN8Mf6a0JdLPzchkvXcfEsFArx0UcfAXDdddcxefJkrrrqKjZu3MhPfvITVq1a5bmwbqFQiCVLlvDiiy/ywgsvsH//fu68806ys7Pp1KmT77w6paamsmrVKuLi4ti5cyd33XUX+fn5/Pvf/2b06NG8//77vhPr1Lx5c2677TZiYmJqLJs/fz779+8/43WPGzeu1sedczz33HOUlJSc8bojIS0tjRUrVnDgwAE6dOjAp59+SlJSErt27eKHP/xh1c9cNNJrgj/WXxOs/9yGS5+4eHbim05RURFXXXUVAF26dOHo0aO+ssIWExNDYmIi48aNY/369bzyyits376djIwM+vXr5zuvTs2bNycuLg6Aiy++mK+//hqAjIyMer3pR0pKSgrjx4/nhRdeqPErISGhXut+6qmniI2NJS4urtqv+Pj4WgelaNOkSRPi4+Np06YNnTp1IikpCaj47xzt/XpN8Mf6a4L1n9twNfMdcLbbtm0b48aNwznH7t27KS8vp2nTpgCUl5d7rju1b39g17dvX/r27cuMGTN49dVXPVWFp3PnzkyaNInBgwczb948MjIyADh69CiHDx/2XHdqv/zlL0/aOXXq1Hqtu2fPntxyyy307NmzxrLZs2fXa92RcOLPzqRJk6oti/b/tnpN8Mf6a4L1n9tw6V5Fnh06dIgWLVrQokULMjIy6N69O3FxcWzbto3169czbNgw34l1Wr9+PTfccEONx1u0aEEoFPJQFL5+/fqRm5vLrFmzuOiii3jyySeJjY1l7969pKWl0blzZ9+JdQqFQiQmJta67PLLL6/Xutu1a0ebNm244IILaizLzMwkOTm5XusPmnOOyy67jJYtW5Kamlr1+IYNG9ixYweDBw/2WFc3vSb4c+JrwoUXXsiTTz7JOeecY+Y1wfrPbbh0jIuIYSUlJRw7dozzzz+fPXv2sGzZMlJTU+nWrZvvNBGRQGhXURQ4fPgwu3fv5jvf+U61x9etW0f37t09VYXPcr/l9ldffZV77rmHJk2a8Nxzz/Hoo4/Stm1b1qxZwzPPPMPw4cPrtf6CggJefvlltm7dCkBSUhIjRowwMxRZ7rfcDrb7LbfXZfHixWRlZfnOaBA6ONezf/7zn7Rp04aUlBQyMzOrXavhxNMio9Xx/tTUVHP91rf91KlT2bBhA++//z7Z2dn89a9/5a233mL58uX1Psbl2WefZcCAARw8eJDMzEwyMjI4ePAgAwYM4Nlnn22gf0FwLPdbbgfb/ZbbT6UxnQ6NE6969+7t1q5d644dO+Zmz57t2rdv7z755BPnnHPp6eme607Ncr/ldueqN3br1u2ky85Ely5dXHFxcY3Hi4qKXOfOneu17kiw3G+53Tnb/ZbbnXPu4YcfrvXX2LFjXevWrX3nNRjtKvKsrKys6kDKkSNHkpycTFZWFgsXLjRx+prlfsvtUP0Mk/vvv7/asvqeNnv8uJlvS0hIMHFBMcv9ltvBdr/ldqg4HXrcuHFVZ6GdyMJrWrg0uHhWVlZGWVkZLVu2BKB///68+OKLDB061MTpd5b7LbcDDB06lH379tG6dWvGjBlT9fiGDRvqfaGvG264geuvv557772X9u3bExMTw+bNm5k1a1atZ4xEG8v9ltvBdr/ldjh7TofWWUWeTZgwge9973sMGDCg2uPvvvsuo0aNoqCgwFNZeCz3W24PmnOOuXPn8uqrr1Y7SPGWW27hjjvuoEmT6D48znK/5Xaw3W+5HSA/P5+uXbvWetr2O++8w7XXXuuhquFpcBFpRMaMGcPTTz/tO0NEJDDRPT6epU782N8iy/2W24FAbwR38803B7buSLDcb7kdbPdbbgf7/bXR4BKFrN+F1HK/5Xaoebn1hrRp06bA1h0Jlvstt4PtfsvtYL+/NhpcopD1vXeW+y23A8ybNy+wdVvfNpb7LbeD7X7L7WC/vzY6xiUKrV+/npSUFN8ZZ8xyv+X2E61du5a0tLQGXWdJSQnnnXdeg64zkiz3W24H2/2W28F+f230iUsUOv7GuXbtWs8lZ8Zyv+X2Ew0ZMqTe6zh8+DBfffVV1ffHX/zWrVtX73VHguV+y+1gu99yO9jvD4cGlyjWEG8+Plnut9A+c+bMWn8988wzHDx4sF7rtnwrB7Ddb7kdbPdbbgf7/eHSBeg8mzlzZq2PO+fq/eYTCZb7LbcDPPDAA9x22221XhGzvhfQGz9+PMuWLaNnz57MmTOH6667jsWLF9OjRw8T+8wt91tuB9v9ltvBfn+4NLh4FuSbTyRY7rfcDhW7tcaPH1/rXWv//ve/12vd1m+HYLnfcjvY7rfcDvb7wxbwvZDkFHr27Ok2bNhQ67J27dpFuOb0We633O6cc3PmzHEff/xxrcteeumleq07JSXFlZaWVnts2bJlLikpybVp06Ze644Ey/2W252z3W+53Tn7/eFqOnHixIm+h6ezWWxsLImJiSQmJtZYlpiYWDU9RyvL/ZbbAUKhUK3tQL3bv/rqK5xz1S4dnpycTGZmJn/729+i/kJ9lvstt4PtfsvtYL8/XDodWsSwkpISFi5cWO2+KkOHDiUhIcFzmYhIMDS4RAHrbz6W+y235+XlMXr0aK699lrat2+Pc44tW7bw7rvvMnPmTG666aZ6rd/ytgHb/ZbbwXa/5Xaw3x8OnQ7tWV5eHpdddhn5+fns37+fffv28eabb5KSkkJeXp7vvFOy3G+5HSrOIFi5ciUvv/wyv//973n88cd55ZVXeO+99xg/fny91m1921jut9wOtvstt4P9/rB5OrZGKnXr1s1t2rSpxuNffPGF69atW+SDTpPlfsvtzjnXuXPnM1oWDuvbxnK/5XbnbPdbbnfOfn+49ImLZ+Xl5SQnJ9d4vEOHDpSXl0c+6DRZ7rfcDtCrVy+ys7P58MMP2b17N0VFRXz44YdkZ2eTmZlZr3Vb3zaW+y23g+1+y+1gvz9cGlw8C/LNJxIs91tuB3j++efp0KEDd955Jx07dqRjx47ceeedtG/fnjlz5tRr3da3jeV+y+1gu99yO9jvD5vvj3zOdocOHXKTJ0923bt3d+eee65r3bq16969u5s4caI7ePCg77xTstxvuT1o1reN5X7L7c7Z7rfc7pz9/nDprCKRRmTMmDE8/fTTvjNERAKjXUVRyPpFgiz3W24HWLFiRWDrtr5tLPdbbgfb/ZbbwX5/bTS4RKEg33wiwXK/5XYg0BupWd82lvstt4PtfsvtYL+/NhpcopD1vXeW+y23A8ybNy+wdVvfNpb7LbeD7X7L7WC/vzY6xiUKrV+/npSUFN8ZZ8xyv7X2SF4l09q2+TbL/ZbbwXa/5Xaw31+bZr4DpPY3n0svvdTMJZot91tur+2S/2+++Sa/+tWvArvkv5VtA7b7LbeD7X7L7WC/PxzaVeSZ9Us0W+633A665H9dLPdbbgfb/ZbbwX5/2Dycgi0nsH6JZsv9ltud0yX/62K533K7c7b7Lbc7Z78/XPrExTPrl2i23G+5HXTJ/7pY7rfcDrb7LbeD/f5waXDxzPolmi33W26H2i/5f9ddd+mS/9jut9wOtvstt4P9/rD5/sjnbFfbJZp79Ohh5hLNlvsttwfN+rax3G+53Tnb/ZbbnbPfHy6dDi3SCO3Zs4fzzz/fd4aISIPTrqIotmfPHt8J9WK530L72rVr6datG+eccw4333wzu3fvrlrWv3//wP5eC9umLpb7LbeD7X7L7WC//0QaXDzz9ebTUCz3W24HeOCBB5g+fTqFhYWkpqZyzTXXsG3bNqD+V8u0vm0s91tuB9v9ltvBfn+4NLh4FuSbTyRY7rfcDrBv3z4GDx7MhRdeyKOPPsqECRP4wQ9+wJdffklMTEy91m1921jut9wOtvstt4P9/rD5ObRGjguFQtW+nzt3ruvatavbunVrjWXRyHK/5XbnKq7ZUF5eXu2xl19+2XXp0sUlJSXVa93Wt43lfsvtztnut9zunP3+cOmS/54dOnSIY8eO0aRJxYdft99+O82bN6d///6UlZV5rjs1y/2W2wH69u1Lfn4+WVlZVY+NGDGCmJgYbr/99nqt2/q2sdxvuR1s91tuB/v9YfM9OZ3tsrOz3aJFi2o8/sorr7jmzZt7KDo9lvsttwfN+rax3G+53Tnb/ZbbnbPfHy6dDi1i2LFjx1i+fHm1G6pdffXVVf/HJSLS2GhXURSw/uZjud9y+4oVK7jtttto06ZN1d2ht2zZwo4dO3jppZe4+uqr67V+y9sGbPdbbgfb/ZbbwX5/OPSJi2dBv/kEzXK/5XaAyy+/nDlz5nDFFVdUe3zVqlVkZ2fzySefnPG6rW8by/2W28F2v+V2sN8fNh/7p+T/9OzZ061atarG4x988IHr0aOHh6LTY7nfcrtzznXp0uWMloXD+rax3G+53Tnb/ZbbnbPfH67G89mRUaWlpTX+jxkqbpZl4Shwy/2W2wE6derE5MmTKSoqqnqsqKiISZMm0aFDh3qt2/q2sdxvuR1s91tuB/v94dLg4lmQbz6RYLnfcjvAX/7yFzZv3kxycjLx8fGce+65JCcns2XLFubOnVuvdVvfNpb7LbeD7X7L7WC/P2y+P/I52+3cudPdfffdLj4+3sXFxbn4+HgXHx/v7r77brdjxw7feadkud9y+7cVFRW5oqLPtOHeAAAEOUlEQVSiBluf9W1jud9yu3O2+y23O2e/P1w6ODeKFBcXA3DBBRd4Ljkzlvsttm/ZsoWcnBw2bdrE0KFD+d3vfkdsbCwAffr0YeXKlQ3y91jcNiey3G+5HWz3W24H+/110a4iz7Zs2cKAAQPo2rUrU6ZMoVWrVlXL+vTp47EsPJb7LbcD3HfffQwdOpTc3Fx27dpF//792b9/P1Cxr7s+rG8by/2W28F2v+V2sN8fLg0ungX55hMJlvsttwNs376d+++/n8zMTF588UUGDx5M//79KSkpqfdNFq1vG8v9ltvBdr/ldrDfHzbf+6rOdt++8dVjjz3mevXq5fbu3WvipliW+y23O1dxk8VvmzZtmsvMzHSdO3eu17qtbxvL/ZbbnbPdb7ndOfv94dKVcz07dOhQte9//etf06JFi2qTcjSz3G+5HSAlJYUlS5YwcODAqsfGjh1LkyZNGDt2bL3WbX3bWO633A62+y23g/3+sPmenM52w4YNc2+99VaNx5988kkXExPjoej0WO633O6cc6Wlpa60tLTWZYWFhfVat/VtY7nfcrtztvsttztnvz9cOqvIs+MXBWrZsmWNZdu2baNt27aRTjotlvsttwfN+rax3G+5HWz3W24H+/3h0uAiIiIiZuisIhERETFDg4uIiIiYocFFRKJCcnIyl112GWlpaXTp0oUbb7yR995775R/7s9//jMFBQURKBSRaKDBRUSixvz581m7di0bN24kOzubQYMG8a9//avOP6PBReTsosFFRKLSjTfeyOjRo3niiSdYunQpffr0IRQK0aNHD1544QUAZs+ezerVq/n5z39Oeno6+fn5ADzxxBN897vfJSMjg0GDBvHll1/6/KeISAPSBehEJGr16tWL119/nYyMDJYvX07Tpk0pLi4mIyODgQMHMmrUKF566SXGjh1LVlYWAPPmzaOgoICVK1fStGlT5s6dy89+9jMWLlzo+V8jIg1Bg4uIRK3jV2soKipi5MiRFBQU0KxZM3bv3s26deu49NJLa/yZ119/ndWrV5OZmQlAeXk5TZs2jWi3iARHg4uIRK1Vq1bRo0cP7r33XoYMGcJrr71GTEwMGRkZJ71pnHOO3/zmN2RnZ0e4VkQiQce4iEhUWrhwIc8++ywPPvgge/bsoX379sTExPDuu++ydu3aque1bt2akpKSqu+HDh3KzJkzKS4uBuDIkSN89NFHEe8XkWDoExcRiRo/+tGPaNmyJQcPHiQ1NZX8/HyuvPJKpk6dyujRo5k6dSqpqan07t276s/k5OTw0EMPMW3aNKZMmcIdd9xBUVER3//+94mJieHo0aOMHDmSUCjk8V8mIg1Fl/wXERERM7SrSERERMzQ4CIiIiJmaHARERERMzS4iIiIiBkaXERERMQMDS4iIiJihgYXERERMUODi4iIiJihwUVERETM+F+yANlGPIh6MgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# graph\n",
    "import pandas\n",
    "from matplotlib.pyplot import figure\n",
    "figure(num=None, figsize=(8, 6), dpi=80, facecolor='w', edgecolor='k')\n",
    "import matplotlib.pyplot as plt\n",
    "from pandas.plotting import register_matplotlib_converters\n",
    "register_matplotlib_converters()\n",
    "plt.plot(df['date'], df['R0'], color = 'magenta',label = 'R0', marker = 's')\n",
    "\n",
    "plt.legend()\n",
    "plt.suptitle('Estimation of R0')\n",
    "plt.title('API from BlankerL and data from Ding Xiang Yuan')\n",
    "plt.grid()\n",
    "plt.xlabel('Date')\n",
    "plt.ylabel('R0')\n",
    "plt.xticks(rotation = 90)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "64269.5 2.7235475907238937\n"
     ]
    }
   ],
   "source": [
    "# sensitivity analysis\n",
    "suspect = 21675\n",
    "confirm = 44762\n",
    "\n",
    "# parameters\n",
    "p = 0.9\n",
    "rho = 0.9\n",
    "Yt = suspect * p + confirm\n",
    "Tg = 10\n",
    "t = 73\n",
    "lamb = math.log(Yt)/t\n",
    "R0 = 1 + lamb * Tg + rho*(1-rho)*(lamb*Tg)**2\n",
    "print(Yt,R0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
